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Abstract 

We investigate structural order in glassy water by performing classical molec- 
ular dynamics simulations using the extended simple point charge (SPC/E) 
model of water. We perform isochoric cooling simulations across the glass 
transition temperature at different cooling rates and densities. We quantify 
structural order by orientational and translational order metrics. Upon cool- 
ing the liquid into the glassy state, both the orientational order parameter 
Q and translational order parameter r increase. At T = K, the glasses 
fall on a line in the Q-t plane or order map. The position of this line de- 
pends only on density and coincides with the location in the order map of 
the inherent structures (IS) sampled upon cooling. We evaluate the energy 
of the IS, eis(T), and find that both order parameters for the IS are propor- 
tional to eis- We also study the structural order during the transformation 
of low-density amorphous ice (LDA) to high-density amorphous ice (HDA) 
upon isothermal compression and are able to identify distinct regions in the 



order map corresponding to these glasses. Comparison of the order parame- 
ters for LDA and HDA with those obtained upon isochoric cooling indicates 
major structural differences between glasses obtained by cooling and glasses 
obtained by compression. These structural differences are only weakly re- 
flected in the pair correlation function. We also characterize the evolution of 
structural order upon isobaric annealing, leading at high pressure to very-high 
density amorphous ice (VHDA). 

I. INTRODUCTION 

Crystalline solids are characterized by a periodic structure and, therefore, by the pres- 
ence of long range-oder. A well-defined theoretical and experimental framework is available 
for characterizing the structure of crystals [1] . On the other hand, numerous materials found 
in nature and widely utilized in technology exist in an amorphous state (e.g., glasses and 
liquids) [2]. Characterizing the structure of amorphous materials is a challenging task be- 
cause such systems organize over short distances and show no long-range order. Moreover, 
the local structure in such materials is not unique and it can change throughout the sample. 
Therefore, new tools are required to characterize the structure of amorphous materials. One 
promising approach is to develop order metrics (order parameters) that quantify certain 
statistical properties of the structure [3]. Recent work in this direction has been performed 
on the hard-sphere [4,5] and Lennard- Jones (LJ) [6] particle systems. Structural order was 
characterized with two metrics: a translational order parameter r quantifying the tendency 
of particle pairs to adopt preferential separations, and a bond-orientational order parame- 
ter Q quantifying correlations between the bond angles defined by a given particle and its 
nearest neighbors. The novelty of this approach resides not in the order metrics themselves, 
which have been used before (e.g., [7]), but in the idea of an order map, in which different 
states are mapped onto a plane whose axes represent the order metrics [4,5]. The order map 
for the Lennard- Jones and hard-sphere systems consists of a line, indicating that the two 
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order metrics are not independent [5,6]. 

Extension of these ideas to liquid water has proven to be most useful. In Ref. [8], the 
structural order of cold liquid water was investigated and found to be closely related to 
the well-known dynamic [9-11] and thermodynamic [12] anomalies of this substance. It 
was found that there exists a region in the phase diagram bounded by loci of maximum 
orientational order and minimum translational order within which the dynamic and ther- 
modynamic anomalies occur. Water's order map was found to be more complex than for 
simpler, spherically-symmetric systems. State points define a two-dimensional region in the 
order map, implying that orientational and translational order are in general independent. 
Remarkably, however, it is found that the region where water exhibits dynamic or thermo- 
dynamic anomalies defines a line in the order map, which implies that when water behaves 
anomalously, its translational and orientational order metrics are no longer independent, 
but are instead strictly correlated. This property was then used to quantify threshold values 
of the order metrics needed in order for various anomalies to occur [8]. Similar (though not 
identical) behavior was found to occur in another tetrahedral liquid, silica [13]. 

Notable among water's peculiarities is the fact that it exhibits polyamorphism [14-16], 
i.e., the presence of more than one glassy form. Uniaxial compression of ice-J^ [17] or 
ice-/ c [18,19] at T = 77 K to pressures P ^ 1 GPa produces a disordered high-density 
material named high-density amorphous ice (HDA). If HDA is recovered at T = 77 K and 
P = 1 bar and then heated isobarically, it transforms irreversibly at T ps 125 K to a dis- 
ordered low-density material named low-density amorphous ice (LDA) [17,20]. Moreover, 
LDA transforms to HDA when compressed at T = 77 K to P = 0.6 GPa [20,21]. The 
LDA-HDA transformation is reversible above T ps 130 K [22]. Furthermore, because it 
is very sharp and shows hysteresis, it has been suggested that the LDA-HDA transforma- 
tion is a first-order transition [22-24]. The LDA-HDA transition plays a central role in 
one of the thermodynamically-consistent scenarios that have been proposed to explain the 
experimentally observed properties of supercooled water [16]. 

There are several other routes to producing amorphous ice [14-16] in addition to those 
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mentioned above. Because glasses are non-equilibrium materials, their properties are history- 
dependent, and the specific procedure followed in the glass preparation may lead to different 
amorphous forms. It is not clear whether all the different amorphous ices can be classified 
into two families [25-27], one corresponding to LDA and the other to HDA. One complication 
standing in the way of such a satisfyingly simple view is the experimental observation that 
HDA, when annealed or heated at high pressure becomes appreciably denser. The resulting 
glass is called very-high density amorphous ice (VHDA) or relaxed HDA (RHDA) [28]. 
Whether VHDA/RHDA is simply yet another form of glassy water or whether, as suggested 
by simulations [29,30], it is the more stable form to which HDA relaxes irreversibly is a 
question currently under investigation. 

It should be evident from the preceding discussion that a precise characterization of 
structural order in water glasses is of considerably interest. The purpose of the present work 
is to provide such a characterization. In particular, (i) we use the order metric approach 
introduced in [8] for equilibrium liquid water to quantify structural changes that occur when 
liquid water is cooled isochorically to the glassy state; (ii) we quantify the changes in struc- 
tural order that accompany the LDA— >HDA transformation; (iii) we investigate structural 
evolution during the annealing/isobaric heating of the amorphous ices obtained upon com- 
pression of LDA; and (iv) we compare the location in the order map of the various glasses 
obtained in (i), (ii) and (iii). By doing this, we are able to characterize quantitatively struc- 
tural differences between samples of amorphous ice obtained by different procedures. In 
recent related work, Guillot and Guissani [31] performed computer simulations of glassy 
water obtained along different paths. They used the radial distribution function and struc- 
ture factor to study the structure of the different glasses. In this work, we instead use 
the order map approach. As we show below, this approach allows a precise and sensitive 
characterization of structure. 

This work is organized as follows. In the next Section, we present the simulations details. 
In Section III, we study the local structure when the liquid is cooled isochorically across 
the glass transition at different cooling rates and for different densities. We study both 
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the structure in the instantaneous configurations and in the corresponding potential energy 
minima attained by steepest descent mapping (inherent structures, IS). In Section IV, we 
study the evolution of structural order during the LDA— >HDA transformation, and upon 
annealing (i.e., isobaric heating) different glasses obtained during the compression procedure. 
The comparison of the local structure of the glasses obtained with the different procedures 
mentioned above is discussed in Section V. We present the conclusions in Section VI. 

II. SIMULATIONS 

We perform classic molecular dynamics (MD) simulations of a (N = 216)-molecule sys- 
tem using the extended simple point charge (SPC/E) model of water [32]. This model has 
been extensively used to study the thermodynamics [10,33] and dynamics [34,35] of liquid 
water and is consistent with experimental facts (e.g., at low temperatures it exhibits density 
anomalies [36] and a diffusivity that increases upon compression [10]). The SPC/E model 
has also been used to study glassy water, and is able to produce glassy states corresponding 
to LDA, HDA and VHDA [30,37]. 

To cool the liquid to the glass state, we perform constant volume simulations at p = 
0.9, 1.0, 1.1, 1.2, 1.3 and 1.4 g/cm 3 . At each density we average the results over 16 
independent simulations (with the exception of p = 1.00 g/cm 3 where we use 32 independent 
simulations). We use periodic boundary conditions and the reaction field method (with a 
cutoff of 0.79 nm) to treat the long range forces. Upon cooling at each density, we change the 
T by rescaling the velocities of the molecules using the Berendsen thermostat [38] . At every 
time step St = 1 fs, we decrease the thermostat temperature by 5T = q c 5t, where q c < 
is the cooling rate. Cooling simulations starting from an equilibrium liquid at T = 300 K 
down to T = K are performed with cooling rates q c = —30, — 10 2 , — 10 3 , — 10 4 and 
— 10 5 K/ns. Similar values for q c have been used in previous simulations at p = 1.00 g/cm 3 to 
study the effect of cooling and heating rates on the glass transition temperature T g = T g (P) 
[39,40]. The glass obtained with a cooling rate q c = —30 K/ns behaves, upon heating at q h = 
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+30 K/ns, as an experimentally slow-cooled glass, and shows no signs of hyperquenching 
effects [40]. On the other hand, the smallest cooling/heating rates accessible nowadays in 
computer simulations are \q\ ~ fO K/ns [31]. 

To study the evolution of structural order during the LDA— >HDA transformation, we 
perform compression simulations at T = 77 and 170 K and average results over 16 inde- 
pendent trajectories at each temperature. We perform MD simulations at constant p for 
intervals of 1 ps and at the end of each interval we increase p by Sp = 5 x 10~ 5 g/cm 3 . These 
p changes are performed by rescaling isotropically the coordinates of the molecular center 
of mass. Therefore, our compression rate is 5p/5t = 5 x 10~ 5 g/cm 3 /ps. This value was 
already used to study the potential energy landscape region sampled during the LDA-HDA 
transformation [37]. We also perform 16 independent compression simulations at T = K; 
at each simulation step we change p by Sp = 5 x 10~ 5 g/cm 3 and then minimize the energy. 
At each step, p is modified by rescaling isotropically the center of mass of each molecule. 

We stress that our results are relative to the cooling/compression rates studied here. 
However, the cooling/compression rates we use are such that the compression run (0.9 < 
p < 1.4 g/cm 3 ) and the slowest cooling run (300 > T > K) both take 10 ns, i.e. both 
numerical 'experiments' are based on the same time scale (as is expected to be the case in 
an experiment). 

The annealing procedure of the glasses obtained upon compression consists of isobaric 
MD simulations where the temperature is increased from T = 77 K (i.e., the compres- 
sion temperature used in the LDA^HDA transformation) up to T = 170 K. We use 
the Berendsen thermostat [38] to fix the T, and at every time step St — 1 fs, we in- 
crease the thermostat temperature by ST = q h St, where qn = +30 K/ns is the heat- 
ing rate. A coupling to an external bath at a given P (analogous to the Berendsen 
thermostat) is used to keep the pressure constant [38]. We perform annealing runs at 
P = -0.55, - 0.17, 0.01, 0.41, 0.68, 0.84, 1.10, 1.38 and 1.90 GPa and average results 
over 16 trajectories. 

We also investigate structural order in mechanically stable configurations sampled during 
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isochoric cooling. These configurations are local potential energy minima called inherent 
structures (IS). To obtain the IS we perform energy minimizations by using the conjugate 
gradient algorithm. We consider the minimization complete when the change in energy 
between two successive minimizations is < 1CT 15 kJ/mol. 

We compute for the glasses the same translational and orientational order parameters r 
and Q introduced in [8] to study structural order in liquid water. The translational order 
parameter r is given by 



where £ = rpl/ 3 is the distance between the oxygen atoms of pairs of molecules, r, divided 
by the mean nearest-neighbor separation at the number density p n = N/V; goo( r ) is the 
oxygen-oxygen radial distribution function (RDF) and £ c = 2.843 is a cut-off distance. In 
an ideal gas, the RDF is equal to 1 and r = 0. In a crystal, there is long-range order and 
9oo( r ) 7^ 1 over l° n g distances, and r is large. 

The orientational order parameter Q is given by 



where ipj^. is the angle formed by the lines joining the oxygen atom of a given molecule and 
those of nearest neighbors j and k (< 4). For the purpose of this calculation, we limit our 
attention to the four oxygen atoms that are closest to a given oxygen atom. This definition 
is a slightly modified version of the metric introduced in Ref. [41]. Equation (2) implies that 
— 3 < Q < 1. In an ideal gas, < Q >= (where < ... > denotes an ensemble average) [8]. In 
a perfect tetrahedral network cos(-?/>j fe ) = —1/3, and Q — 1. Thus, Q measures the degree of 
tetrahedrality in the distribution of the four nearest oxygen atoms around a central oxygen 
atom. 




(1) 




(2) 
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III. STRUCTURAL ORDER UPON ISOCHORIC COOLING 



A. Instantaneous Configurations 

Figure 1(a) shows the evolution of the order parameters Q and r during cooling of an 
initially equilibrated liquid at constant density p = 1.00 g/cm 3 from T = 300 K down to 
T = K. We show the results for different cooling rates q c and include the line obtained in [8] 
separating the accessible and the inaccessible regions for equilibrium liquids. No equilibrium 
liquid states of the SPC/E model exist in the range 220 K< T < 400 K and 0.85 g/cm 3 < 
p < 1.3 g/cm 3 with order parameter combinations (Q,r) to the right of this limiting line. 

The starting location of the system in the ordering map is indicated by the square 
symbol 'A' and is in agreement with the values reported in [8] for a liquid at T = 300 K and 
p = 1.0 g/cm 3 . Upon cooling at the slower rate q c = —30 K/ns, the system moves in the 
order map along the line delimiting the accessible region for the liquid state and finally at 
approximately T = 240 K, the system departs from this line and moves into the accessible 
region. The boundary between the accessible and inaccessible regions is the locus of state 
points where structural, dynamic, and thermodynamic anomalies occur [8]. Equilibrium 
liquid states that display structural, dynamic, or thermodynamic anomalies (decrease of both 
r and Q upon isothermal compression; increase in diffusivity upon isothermal compression; 
decrease in density upon isobaric cooling, respectively) are therefore characterized by only 
one order parameter (line in the Q-t plane). In contrast, liquid states that do not exhibit 
anomalies are characterized by two order parameters (accessible area in the Q-t plane) [8] . 
Therefore, at p = 1.0 g/cm 3 and 300 K > T > 240 K, the equilibrium liquid experiences a 
cascade of anomalies (first structural, then dynamic, finally thermodynamic) upon cooling 
[8]. For T < 240 K, the quenched liquid falls out of equilibrium and in so doing moves into 
the accessible region of the order map. We use the term equilibrium to denote states where 
the structural relaxation time is short compared to the simulated time at the given T and 
p. Equilibrium states for the liquid at p = 1.0 g/cm 3 and 220 K< T < 300 K lie on the 
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limiting line separating the accessible and inaccessible regions [8]. 

In order to understand the trajectory followed by the quenched system at p = 1.00 g/cm 3 
in the order map we note that the liquid is able to reach equilibrium during cooling down 
toT« 240 K. The structural relaxation time at T > 240 K is t r < 30 ps [10] while during 
the cooling process, the typical time scale for the system to relax before T changes by 1 K 
is dt — 1 K/|g c | ~ 33 ps ~ t r . For T < 240 K, therefore, the system is not able to reach 
equilibrium upon cooling (i.e., dt < t r ) and it attains, for a given Q, larger values of r 
than those possessed by the liquid in equilibrium at the given T and p. Equivalently, the 
quenched system attains smaller values of Q, for a given r, than those of the equilibrium 
liquid at the given T and p. 

It is interesting to study the behavior of the thermodynamic properties when the system 
falls out of equilibrium upon cooling. Figure 1(b) shows the behavior of pressure P as a 
function of T upon cooling at q c = —30 K/ns and q c = — 10 4 K/ns. As a guide to the eye, 
we also include data from equilibrium simulations taken from Ref. [33]. Upon cooling at 
q c = —30 K/ns the pressure clearly follows the equilibrium P(T) trend, showing a minimum 
at T ps 250 K. A minimum in P(T) at constant density corresponds to a maximum in p(T) 
at constant density (see e.g. [10]). Therefore, upon cooling at q c = —30 K/ns, the liquid is 
able to experience a density anomaly [(dP/dT) < 0] in the approximate range 250 K> T > 
190 K, whereas for the equilibrium liquid, density anomalies occur for all T < 250 K. When 
the liquid is cooled faster, at q c = — 10 4 K/ns, the pressure deviates from the equilibrium 
P(T) curve at a comparatively higher T. No signature of a density anomaly is observed when 
cooling at this high rate. This suggests that configurations responsible for producing density 
anomalies possess comparatively larger equilibration times. This observation warrants a 
systematic investigation. We note that pressure fluctuations make it difficult to estimate 
the T at which the instantaneous P deviated from the equilibrium P(T) curve. 

The phenomenology discussed above when relating the behavior of the order parameters 
upon cooling at q c = —30 K/ns is common to all the cooling rates studied [see Fig. 1(a)]. 
However, we observe that the larger \q c \, the higher the T at which the system falls out of 
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equilibrium (i.e., the earlier the system deviates from the line delimiting the inaccessible 
region). Furthermore, the larger is \q c \, the closer are the final values of Q and r to those 
of the starting liquid (square symbol 'A' in the figure), meaning that the structure of the 
glass is more similar to that of the starting liquid. Equivalently, the smaller is \q c \, the more 
structured is the final glass, i.e., the larger the final values of Q and r. 

The long-dashed line in Fig. 1(a) is the locus of Q-r values for p = 1.0 g/cm 3 glasses 
obtained at T = K with the five q c studied. Interestingly, such isochorically-quenched 
glasses define a line in the order map. We note that these results are similar to simulations 
performed using a Lenard- Jones potential [43]. In fact, Fig. 7 in that work shows that upon 
isochorically cooling a liquid with different cooling rates, the orientational and translational 
parameters increase monotonically with decreasing cooling rate, and the glasses obtained at 
the lowest T fall on a straight line with positive slope in the order map, just as we find for 
SPC/E water. 

The features observed at p = 1.00 g/cm 3 are also found at all the other densities studied, 
i.e., p = 0.9, 1.1, 1.2, 1.3, and 1.4 g/cm 3 . The results for p = 1.3 g/cm 3 are shown in 
Fig. 1(c). In agreement with [8], comparison between Figs. 1(a) and (c) indicates that the 
location in the order map of the starting liquid at T = 300 K (square symbol) depends on 
the density. Furthermore, at the high density p = 1.3 g/cm 3 the equilibrium liquid does 
not exhibit structural, dynamic, or thermodynamic anomalies. Accordingly, the quench 
trajectories in the order map do not attain the line separating the accessible and inaccessible 
regions. We note that at all the densities and cooling rates studied, the glasses always 
remain out of the inaccessible region obtained previously for the equilibrium liquids [8]. In 
this respect, the inaccessible region in the order map seems to be the result of topological 
constrains inherent in the molecular interactions and does not depend on whether the system 
is in the liquid or glassy state. 

The radial distribution function is one of the standard tools used in experiments, theory, 
and simulations to characterize the structure of condensed matter. We show in Fig. 2 the 
oxygen-oxygen radial distribution function gooij') for the glasses obtained at T = K upon 
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cooling at p = 1.00 g/cm 3 with the same five cooling rates shown in Fig. 1(a). These 
glasses are located in the order map along the long-dashed line. While differences in the 
structures of the glasses are suggested by both Fig. 1(a) and Fig. 2, the differences observed 
in gooij') are quite subtle and could be difficult to distinguish from typical experimental 
error bars. It therefore appears that the location of the system in the (Q, r) order map 
is a sensitive indicator of glass structure. Same conclusion holds when comparing gooij') 
for the glasses obtained upon compression of LDA at T = K (where no annealing is 
possible) and the isochorically fast-quenched glasses at the same T. For example, inset of 
Fig. 2 compares the gooij') for the glass obtained upon isothermal compression at T = K to 
p = 1.0 g/cm 3 and for the isochorically quenched glass at q c = — 10 5 K/ns at the same density 
and temperature. The glasses obtained by isothermal compression and those prepared by 
other routes will be discussed and compared in detail in the rest of the paper. The point 
here, though, is to emphasize that the differences in gooij') shown in the inset of Fig. 2 are 
quite subtle, whereas the two glasses have very different order metrics: Q = 0.83, r = 0.57 
for T = K-compression of LDA; Q = 75, r = 50 for the isochorically-quenched glasses at 
q c = -10 5 K/ns. 



B. Inherent Structures 

We now investigate the location in the order map of the IS sampled by the system upon 
cooling from T = 300 K down to T = K. We denote the tetrahedral and translational 
order parameter of the IS by Qis and tjs, respectively. 

Figure 3(a) shows the evolution of Q and r upon cooling the liquid at p = 1.00 g/cm 3 
at q c = —30 K/ns [taken from Fig. 1(a)]. For each configuration obtained from the MD 
simulations, we find the corresponding IS and then calculate the order parameters Qis and 
t is . The evolution of Qis and r IS upon cooling at q c = —30 K/ns are also indicated in 
Fig. 3(a). Interestingly, we find that the Q IS and t is values fall on the same line that 
characterizes the glasses obtained upon cooling isochorically to T = K at different q c 



11 



[long-dashed line in Fig. 1(a); included also in Fig. 3(a)]. In fact, we find that for all the IS 
obtained at p = 1.00 g/cm 3 from systems quenched at any quenching rate q c , the parameters 
Qis and tis fall on the same line that characterizes the isochorically-quenched glasses at 
T ps K [see Fig. 3(b)]. We also note that the values of the order parameters are larger in 
the IS than in the liquid (i.e., the IS are more structured than the liquid configurations). 

In Fig. 3(c) we show that the behavior found at p = 1.00 g/cm 3 occurs at all densities 
studied. For each density, there is a line in the order map that characterizes the location of 
the IS sampled upon cooling, which coincides with the location of the glasses at T = K. 
The location of each of these lines depends non-monotonically on p. At low p, the IS-line is 
located at high Q and r. As density increases, this line moves toward lower values of Q and 
r, and finally, at high p, it shifts to lower values of Q but larger values of r. The range of 
order parameters sampled by the IS shrinks appreciably upon compression within the range 
of q c values studied. 

From Ref. [8], it is found that equilibrium liquid states of a given density fall on a single 
curve in the (Q, r) plane. Therefore, it is interesting to investigate how this equilibrium 
liquid line compares to the corresponding IS line at the same density. This is done in 
Fig. 3(d) where we show the location in the order map of the equilibrium liquid at different 
p (from Ref. [8]) and the lines shown in Fig. 1(c) for the corresponding p. Clearly, for a 
given p the equilibrium liquid line differs from the corresponding IS line. 

The finding that the Qis- t is values fall on a density-dependent line in the order map has 
the following implications: (i) Qjs and Tjg are not independent for a fixed density, i.e., one 
cannot change the average orientational order in the IS without changing the translational 
order; (ii) the function tis(Qis) is independent of the cooling rate, i.e., on the glass history; 
and (iii) the function tis{Qis) is independent of whether the IS corresponds to a system in 
equilibrium or out of equilibrium. Conclusions (ii) and (iii) are unusual since expressions 
relating properties of the IS sampled by the system in equilibrium do not necessarily apply 
to IS sampled by the system out of equilibrium. For example, the dependence of basin 
curvature on IS depth has been found to be different in IS sampled by equilibrium liquids 
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than in IS obtained from fast-quenched or fast-compressed glasses [40,44], and to depend on 
the cooling/compression rates. 

In the potential energy landscape (PEL) approach, any property of an equilibrium system 
[45] or of a system not too far from equilibrium [40,44,46] is known once the energy (eis), 
local curvature (Sis), and pressure of the IS (Pis) are determined [45]. We investigate the 
connection between these properties and the structural order parameters tis-Qis- To do 
this, Fig. 4(a) shows 1 — Qis as a function of T at p = 1.00 g/cm 3 for all q c (similar results are 
obtained for r/s). This figure is reminiscent of the behavior of eis observed in simulations 
of a Lenard-Jonnes system [47] and in water [40] upon cooling an equilibrium liquid, and 
suggests that 1 — Qis ~ eis (similarly, we find that 1 — T/s ~ e.is) meaning that both Qi S 
and Tis are a measure of the basin depth ei S in the PEL. This is confirmed by Fig. 4(b) 
where we show 1 — Qis as a function of eis- 

IV. STRUCTURAL ORDER OF LOW-DENSITY, HIGH-DENSITY AND 

ANNEALED AMORPHOUS ICE 

A. Low-density and High-density Amorphous Ice 

Next, we study the structural order parameters in glassy water. Specifically, we focus 
on the values of Q and r corresponding to the glasses obtained during the LDA^HDA 
transformation. We quantify the tetrahedrality and translational order of both groups of 
glasses and identify regions in the order map corresponding to LDA and HDA. So far, such 
a characterization has been missing. 

Figure 5(a) shows the p-dependence of pressure P during the LDA^HDA transformation 
at the compression temperatures of T = 0, 77, and 170 K, which lie below the glass transi- 
tion locus T g (P) for this system, as determined by the behavior of C P (T) upon heating (not 
shown). In experiments, LDA is obtained by isobaric heating of HDA at P = 1 atm. In com- 
puter simulations LDA has been obtained by isochoric cooling of a low-density liquid [42,48]. 
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Following [42,48], we obtain the starting LDA by isochoric fast-cooling of an equilibrium 
liquid at T = 220 K. The procedure followed to obtain LDA should not alter qualitatively 
our results (see Appendix). In Fig. 5(a), LDA corresponds to the linear behavior of P(p) 
at low-p. The transformation of LDA to HDA starts at the end of this linear regime and is 
indicated on the left side of the figure. It is not so clear when the transformation to HDA is 
complete, in particular because there are aging/annealing effects [26]. For simplicity, in the 
following analysis we assume that HDA corresponds to the glasses obtained by compression 
whose density exceeds p = 1-3 g/cm 3 . The location of these glasses in the P — p plane is 
indicated in Fig. 5(a). 

Figure 5(b) shows the evolution in the order map of the glasses obtained during the 
LDA^HDA transformation. The isotherms in Fig. 5(b) correspond to those shown in 
Fig. 5(a). The symbols in each case denote the densities 0.9, 1.0, 1.1, 1.2, 1.3, and 
1.4 g/cm 3 . At the three temperatures studied, the orientational order decreases monotoni- 
cally upon compression, i.e., the glasses become less tetrahedral. However, the translational 
order parameter behaves non-monotonically, exhibiting a minimum in the range from 1.25 to 
1.3 g/cm 3 . An increase in translational order upon compression is the expected behavior for 
a condensed-phase system [4-6] . In liquid water, this behavior occurs only at high densities 
[8]. It can be seen from Fig. 5(b) that the same is true in glassy water. 

It is interesting that LDA and HDA are in well-defined and non-overlapping regions of 
the order map. The difference in tetrahedrality is especially notable, with LDA-like glasses 
having Q values > 0.8 and HDA-like glasses having Q values approximately < 0.7. We also 
note that the compressed glasses remain in the accessible region of the order map. 

B. Annealed Amorphous Ice 

Recently, it has been found that upon annealing or heating isobarically samples of HDA 
at high-P, the glass becomes denser [49]; the resulting glass is called very-high-density 
glass (VHDA) or relaxed HDA (RHDA) [26,29,30]. VHDA/RHDA can also be obtained in 



14 



computer simulations using the SPC/E [30] and the TIP4P models [29]. In fact, it has been 
found [26] that annealing effects are present upon heating not only HDA, but also LDA 
and intermediate glasses in the LDA^HDA transformation. Figure 6(a) shows P(p) during 
compression of LDA at T = 77 K and T = 170 K. We also indicate the effect of isobarically 
heating nine different glasses from T = 77 K up to T = 170 K [26] at pressures ranging from 
P = —0.55 up to 1.9 GPa (in the case of P = —0.55, the maximum T reached is T = 160 K). 
Upon heating, the density of the compressed glasses shifts from the (T = 77 K)-isotherm 
(open circles) to the full diamonds located near the (T — 170 K)-isotherm. 

The changes in structural order that occur during annealing/isobaric heating of the nine 
glasses indicated in Fig. 6(a) are shown in Fig. 6(b). The effect of annealing LDA is mainly 
to reduce the translational order while only slightly reducing the orientational order. This 
is clear from the isobaric heating of LDA, at P = —0.55, — 0.17, and 0.01 GPa [trajectories 
starting from the first three circles (high-Q) of the 77 K-isotherm in Fig. 6(b)]. However, 
as P is increased approaching the LDA^HDA transformation region, the tetrahedral order 
becomes more sensitive to annealing. This is clear from the heating at P = 0.41 GPa, which 
corresponds to a glass located at the end of the LDA region in Fig. 6(a) [trajectory starting 
from the fourth circle in Fig. 6(b)]. Although the order parameters indicate a considerable 
change in structure, in the case of the isobaric heating at P = 0.01 GPa, the density of the 
glass barely changes. 

In contrast to the behavior observed upon isobaric heating of LDA, annealing HDA 
mainly reduces the orientational order, while small changes in r are found to occur. This 
is clear from the isobaric heatings at P = 1.38 and 1.9 GPa. [trajectories starting from 
the last two circles (low-Q) of the 77 K-isotherm in Fig. 6(b)]. We note that in the case of 
the isobaric heating at P = 1.9 GPa, there is actually an increment in translational order, 
the opposite effect to that observed upon heating LDA. In the case of isobaric heating at 
intermediate pressures (approximately 0.41 GPa < P < 1.38 GPa), both the translational 
and orientational order parameters decrease appreciably. 

We stress that upon annealing HDA (P = 1.38 and 1.9 GPa), i.e. during the 
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HDA^VHDA transformation, the location of the glass in the order map shifts to a 
\ow-Q subregion already included in the HDA domain. Similarly, annealing LDA (P = 
—0.55, — 0.17, and 0.01 GPa) shifts the location of the glass in the order map to a low-r 
subregion already included in the LDA domain. This does not apply to glasses of interme- 
diate densities (0.98 g/cm 3 < p < 1.25 g/cm 3 , approximately) which cannot be character- 
ized as LDA or HDA before annealing. These systems start from, and evolve toward, an 
intermediate-Q region in the order map 0.7 < Q < 0.8). 

A salient feature of Fig. 6(a) and 6(b) is the similarity, both in density and in structural 
order, between isobarically annealed glasses at T = 170 K and isothermally (T = 170 K) 
compressed glasses at the same density. 

V. COMPARISON OF STRUCTURAL ORDER IN ISOTHERMALLY 
COMPRESSED, ISOBARICALLY ANNEALED, AND ISOCHORICALLY 

COOLED GLASSES 

A. Isochorically Cooled Glasses and the Low-density and High-density Amorphous 

Ices 

We first compare structural order in LDA and HDA with that in isochorically-quenched 
glasses (Section III) at the same density. To avoid aging effects, we compare the glasses 
obtained upon compression of LDA at T = K to the isochorically-cooled glasses obtained 
also at T = K. The location of these isochorically-quenched glasses in the P — p plane is 
indicated by open symbols in Fig. 7(a). Full symbols denote the glasses from the LDA^HDA 
isothermal compression at T = K. Clearly, glasses obtained upon compression have very 
different pressures than isochorically-cooled glasses of the same density over the range of 
compression and cooling rates explored here. In fact, we find that the glasses obtained by 
compression of LDA (at the present compression rate) at a given p are less stable than those 
obtained by slow cooling at the same density. In other words, e/s for the compressed glasses 
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is higher than for the isochorically cooled glasses at q c = —30 K/ns. 

These glasses are also different structurally. In Fig. 7(b) we compare the evolution of the 
isothermal LDA^HDA compression at T = K in the (Q, r) plane with the corresponding 
location of isochorically quenched glasses. The LDA— >HDA compression is represented by 
the irregular curve. For the isochorically-quenched glasses we show the (Q, r) coordinates 
of the slowest- and fastest-cooled glasses (q c = —30 K/ns, high-Q end; q c = — 10 5 K/ns, 
low-Q end, respectively), joined by a straight line, along which lie the glasses cooled at 
intermediate rates [see, e.g., Fig. 1(a)]. The symbols denote different densities, as indicated 
in the figure. 

Among the low-p isochorically-cooled glasses (p = 0.9 and p = 1.0 g/cm 3 ), only those 
cooled at the slowest rates (q c = —30 and q c = — 10 2 K/ns) fall inside the LDA region. For 
\q c \ > 10 2 K/ns we find that the structure corresponds to glasses lying between LDA and 
HDA in the (Q,r) plane. As previously shown [40], glasses obtained at rates \q c \ > 10 2 K/ns 
are highly metastable and as soon as they are heated, they relax toward those glasses 
obtained by slower cooling rates. 

At high densities, p > 1.2 g/cm 3 , isochorically-quenched glasses do not fall into the HDA 
region in the order map, as one would expect based on their density. Even the slow-cooled 
glasses are located in the region between that of LDA and HDA. This again highlights 
the fact that the glasses obtained by compression in the LDA— >HDA transformation are 
structurally quite distinct from those obtained by isochoric cooling (at least with the cooling 
and compression rates accessible with the present computer times). 

In fact, we have found (results are not shown here) that e/s, Sis, an d Pis of the glasses 
obtained upon isochoric cooling at a given density and those obtained upon compression 
of LDA at the same density can be very different and are very sensitive to the cooling 
rate. A similar situation seems to happen when comparing glasses obtained upon isobaric 
cooling with those from the LDA— >HDA transformation. For example, VHDA/RHDA can 
be obtained upon isobaric cooling only upon fast-cooling while HDA is not accessible upon 
isobaric cooling of the liquid [30]. 
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B. Isochorically Cooled Glasses and the Annealed Amorphous Ices 



Figure 8(a) shows the location in the P-p plane of the annealed/isobarically heated 
glasses indicated in Fig. 6(a) (diamond symbols) and of the glasses cooled isochorically down 
to T = 170 K at different cooling rates (open symbols). Comparison of Fig. 7(a) and 8(a) 
shows that at this relatively high temperature, equation-of-state differences between glasses 
formed by different routes, though still evident, are much less pronounced. 

The location in the order map of the glasses indicated in Fig. 8(a) is shown in Fig. 8(b). 
For comparison, we also include the location of the glasses obtained by isothermal compres- 
sion at T = 170 K (long-dashed line). The annealed glasses are structurally very similar 
to those obtained upon isothermal compression, i.e., these glasses fall on the T = 170 re- 
compression curve in Fig. 8(b). 

We also show in the same figure the "isotherms" connecting glasses formed by isochoric 
cooling at the same rate. Symbols identify the quench rate, and along each "isotherm" the 
density increases monotonically from 0.9 g/cm 3 (high-Q end) to 1.4 g/cm 3 (low-Q end) in 
0.1 g/cm 3 increments. Note that neighboring points along an "isotherm" are not connected 
by a continuity of states, but were instead attained by isochoric quenches at different den- 
sities. In all these cases, Q decreases monotonically upon compression, while r decreases 
initially up to p ~ 1.0 g/cm 3 and increases thereafter. With exception of the isotherms cor- 
responding to the fastest cooling rates (q c = — 10 4 , — 10 5 K/ns), the low-density portions 
(p < 1.25 g/cm 3 ) of all the curves collapse onto a common isotherm. We thus find that 
glasses with 0.9 g/cm 3 < p < 1.3 g/cm 3 , when sufficiently hot, share with the cold equilib- 
rium liquid the remarkable feature [8] of having only one independent order parameter, it 
being impossible to change one (e.g., Q) without changing the other. 

We note, finally, that at high enough densities, appreciable structural differences between 
glasses formed by different routes are evident even when their thermal properties are very 
similar. For example, all the isochorically cooled glasses at p = 1.3 g/cm 3 are at almost the 
same P [see Fig. 8(a)]. However, Fig. 8(b) shows that these glasses have very different order 
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parameters (next-to-last open symbols on each isotherm, on the low-Q end). 

VI. CONCLUSIONS 

In this work we have used order metrics Q and r and the order map concept [4-6,8], to 
quantitatively characterize structural order in glassy states of the SPC/E model of water. 
Our main findings are recapitulated below. 

At a given density, glasses formed by isochoric cooling down to T = K define a line in 
the order map. This line shifts non-monotonically across the order map upon changing the 
density, and corresponds also to the location in the (Q, r) plane of the IS sampled by the 
system both in and out of equilibrium. The finding that the order parameters of the IS fall 
(for a given density) on a line in the order map implies that Q and r are not independent. 
In fact, we find that both Q and r are linear functions of the depth of the IS, i.e. ejs- 

We have also investigated in detail the structural order of the glasses formed upon isother- 
mal compression across the LDA^HDA transformation. At the three temperatures studied 
(0, 77 and 170 K), the orientational order decreases monotonically upon compression. How- 
ever, the translational order decreases until p ~ 1.3 g/cm 3 (approximately, the HDA forma- 
tion density); further compression increases the translational order. We were able to identify 
two non-overlaping and well-defined regions: a high-Q and high-r region corresponding to 
LDA, and a low-Q and low- or intermediate-r region, corresponding to HDA. 

Annealing (i.e., isobaric heating) of compressed glasses is important experimentally: it 
has been found that it can change dramatically their properties [25,49,50]. Accordingly, we 
also have studied the evolution of structural order in isobarically-annealed glasses over a 
wide range of P. We find that annealing LDA not only decreases the density but also causes 
a significant decrease in translational order, while the orientational order barely changes. 
Annealing HDA (this corresponds to the HDA^VHDA/RHDA transformation [29,30,49]), 
reduces the orientational order appreciably, and causes a slight increase in translational 
order. While structural changes upon annealing are significant, they do not cause LDA or 
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HDA to move outside of their respective regions in the order map (defined before annealing). 

We compared the structure of glasses formed by different routes in order to test whether 
we can relate the LDA/HDA glasses with those obtained by isochoric cooling at the same 
density (at least, with the present cooling/compression rates available in computer simula- 
tions). The structure of the glasses formed by different routes can be compared by traditional 
methods, e.g., the oxygen-oxygen radial distribution function goo( r ), or, alternatively, by 
the order map characterization utilized here. When comparing the structure of isochorically 
cooled, and compressed glasses at same T and p, we find that the location in the order map 
is very sensitive to structural changes. Among isobarically-cooled glasses at T = K and 
p = 0.9 or 1.0 g/cm 3 , only those cooled at the slowest rates (q c = —30 and — 10 2 K/ns) 
fall into the LDA region in the order map. In the case of the isochorically-cooled glasses 
at T = K and p = 1.3 or 1.4 g/cm 3 , they fall outside of the HDA region in the order 
map. Therefore, isochorically-cooled glasses in the appropriate density range cannot always 
be classified as LDA or HDA. 

In contrast to the pronounced path-dependence of the structure of 'cold' glasses, we find 
much weaker history dependence in the properties of glasses formed after isobaric annealing 
at T = 170 K, when compared with those attained by isochoric cooling glasses down to the 
same T. Both the annealed and slow-cooled glasses are structurally very similar, and they 
are located in the order map on the same isotherm as compressed glasses at T = 170 K. 
We also find that glasses formed by isochoric cooling to T = 170 K collapse onto a common 
'master' isotherm in the order map regardless of their density as long as \q c \ < 10 4 K/ns. 

As a final observation, we note that all the glasses obtained in this work, under any 
condition and procedure fall in the accessible region of the order map, as found in [8] for 
liquid water. This suggests that the inaccessible region (at least for systems lacking long- 
range order) is a consequence of constrains inherent to the molecular interactions and is not 
dependent on whether the system is a liquid or a glass. 

In this work we have applied the order map methodology to the investigation of structural 
order in water glasses. We find that location in the order map is a sensitive descriptor of 
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structure. We further find that at low enough temperature, glasses obtained by isochoric 
cooling are very different from those formed by isothermal compression, even when compared 
at the same density and temperature. We were able to identify distinct regions of the order 
map corresponding to LDA (Q > 0.8, r > 0.46), and HDA (Q < 0.7, 0.43 < r < 0.5), and to 
follow the evolution of structural order upon isobaric annealing. These results demonstrate 
the usefulness of an order metric approach in the quest for a more precise characterization 
of structural order in amorphous systems in general and water glasses in particular. 
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APPENDIX 

To test how sensitive the simulation results are to the method of preparation of LDA, 
we compare the compression of LDA prepared by isochoric cooling of the equilibrium liquid 
at p = 0.9 g/cm 3 , with that of LDA prepared by isobaric cooling of equilibrium liquid 
at P = 0.1 GPa. In fact, this last procedure is the recipe followed in experiments to 
obtain hyperquenched-glassy water, HGW [51]. X-ray and neutron diffraction measurements 
suggest that LDA is structurally identical to HGW [52] and, although small differences have 
been found [53-55], the common view at present is that HGW and LDA are the same 
material [14,16,56]. 

Figure 9(a) shows the p-dependence of P when compressing both samples of LDA at 
T = 77 K. Both P(p) curves behave in a qualitatively similar way. It can be seen that the 
two compression curves are similar, especially at high enough densities. Equation-of-state 
properties reveal some sensitivity to sample history at low densities. 
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Figure 9(b) shows the corresponding trace in the order map of the two compressions. We 
observe that both samples of LDA exhibit a very similar structural evolution. Therefore, 
Figure 9 suggests that our results do not depend strongly on the way in which LDA is 
prepared. 
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FIGURES 
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FIG. 1. (a) Behavior of the orientational and translational order parameters (Q and r, re- 
spectively) when cooling an initially equilibrium liquid (point 'A': T = 300 K, p = 1.0 g/cm 3 ) at 
constant density p = 1.00 g/cm 3 . Upon cooling at the slower rate, q c = —30 K/ns, the system 
initially moves along the path corresponding to equilibrium Q — r values (in (a) , this path coincides 
with the dot-dashed line delimiting the inaccessible region for the liquid state). At lower T, the 
system deviates from the equilibrium path. The long-dashed line is the result of interpolating the 
Q-t values of the glasses obtained at T = K at different q c . The various cooling rates correspond 
to q c values of -10 5 , - 10 3 , - 10 3 , - 10 2 and -30 K/ns (left to right), (b) Pressure as a function 
of T upon cooling the liquid at p = 1.00 g/cm 3 for two different cooling rates, (c) Same as (a) for 
p = 1.30 g/cm 3 . 
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FIG. 2. Radial distribution function goo( r ) corresponding to the glasses at T = K prepared 
by isochoric cooling at p = 1.00 g/cm 3 and different cooling rates. As the cooling rate increases 
(see arrows), mild changes in the first minimum and the second maximum of goo( r ) are observed 
indicating that molecules in the second shell move closer to the central molecule. These glasses show 
a very similar goo(i") but have quite different order parameters [see Fig. 1(a)]. Same conclusion 
holds when comparing goo( r ) for the glass obtained upon T = K-isothermal compression at 
p = 1.0 g/cm 3 and for the isochorically quenched glass at q c = — 10 5 K/ns at the same density 
and temperature (see inset). The cutoff value £ c = r c pl/ 3 used in the definition of r (see Eq. 1) 
corresponds to a cutoff distance r c = 0.883 nm at p = 1.0 g/cm 3 . The value of r c is indicated by 
the arrow in the lower right corner. 
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FIG. 3. (a) Order parameters Q and t for the slowest cooling rate shown in Fig. 1(a) (square 
symbols). For each configuration obtained upon cooling (e.g., points A, B, and C), we obtain the 
corresponding IS (points A', B', and C). The order parameters in the IS (Qis and r/5), coincide 
with the long-dashed line obtained in Fig. 1(a) for the glasses quenched at T = K. (b) Location 
in the order map of the IS sampled upon cooling at different cooling rates. All IS fall on the 
same long-dashed line shown in (a), (c) Comparison of the location in the order map of the IS 
sampled at different densities at the slowest cooling rate studied, q c = —30 K/ns. There is a line 
characterizing the location of the IS in the order map at each p. This line shifts in the order map 
non-monotonically with p. (d) Comparison between the IS lines shown in (c) and the lines obtained 
by the equilibrium liquid state points at different p [8] (shown at the lower left corner). At each 
density, the shape of the symbols (e.g., left-pointing triangles at p = 1.3 g/cm 3 ) is the same for 
the IS and the equilibrium liquid. 
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FIG. 4. (a) Orientational order parameter 1 — Qj$ in the IS as a function of temperature upon 
cooling at p = 1.00 g/cm 3 at different cooling rates. 1 — Qjs(T) behaves similarly to the energy of 
the IS observed in [40] under same conditions, (b) Qjs as a function of the energy of the IS, ejg, 
suggesting that 1 — Qjs ~ ejs, independently of the cooling rate. Similar results are found at all 
densities studied. 
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FIG. 5. (a) Evolution of pressure with density upon compression of low-density amorphous ice 
(LDA) to produce high-density amorphous ice (HDA). Compression is performed at T = 0, 77, 
and 170 K, below the glass transition temperature. We indicate approximately the location cor- 
responding to LDA and HDA. (b) Order parameters of the glasses sampled along the isotherms 
shown in (a). The location in the order map corresponding to LDA and HDA is also indicated. 
Filled symbols correspond to the glasses obtained at p = 0.9 up to 1.4 g/cm 3 in steps of 0.1 g/cm 3 . 
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FIG. 6. (a) Evolution of pressure with density upon compression of low-density amorphous 
ice (LDA) to produce high-density amorphous ice (HDA). Compression is performed at T = 77, 
and 170 K. Selected glasses at T = 77 K (open circles) are annealed/isobarically heated up to 
T ~ 170 K (filled diamonds); see horizontal arrows. The annealed glasses in the P — p plane 
approach or cross the T ~ 170 K-isotherm. (b) Location in the order map of the glasses indicated 
in (a). Density increases from right to left along both isotherms. 
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FIG. 7. (a) Comparison of the location in the P — p plane of the glasses at T = K obtained 
by isothermal compression of LDA and by isochorically cooling at different densities and cooling 
rates. At T = K no annealing is possible. Open symbols represent different cooling rates while 
solid symbols on the T = K-isotherm correspond to glasses from p = 0.9 g/cm 3 to p = 1.4 g/cm 3 
in steps of 0.1 g/cm 3 . (b) Location in the order map of the glasses indicated in (a). Glasses formed 
by isochoric cooling lie on the straight lines connecting fast- (low-Q end) and slow-cooled (high-Q 
end) extremes at each density. We also show the location of the glasses along the T = K-isotherm 
indicated in (a). 
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FIG. 8. (a) Comparison of the location in the P— p plane of the glasses at T = 170 K obtained by 
isochoric cooling at different densities and cooling rates (open circles), and by annealing/isobarically 
heating of glasses obtained by compression of LDA at T = 77 K (solid diamonds) [see Fig. 6(a)]. 
(b) Location in the order map of the glasses indicated in (a). Lines connecting open symbols 
correspond to glasses obtained at a given cooling rate for different densities [density increases from 
0.9 g/cm 3 (high-Q values) to 1.4 g/cm 3 (low-Q values), in steps of 0.1 g/cm 3 ]. For comparison, 
we also show the location of the glasses obtained by isothermal compression (T = 170 K) of LDA 
(long-dashed line). We note that all glasses in (b) are at T = 170 K and, therefore, they have the 
same vibrational contributions. 
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FIG. 9. Comparison of the /^-dependence upon compression for two samples of LDA prepared 
in a different way: LDA obtained by isochoric fast-cooling of equilibrium liquid at T = 220 K and 
p = 0.9 g/cm 3 , and LDA obtained by isobaric slow-cooling at P = 0.01 GPa. (b) Comparison of 
the order parameters in the glasses sampled upon compression both samples of LDA shown in (a). 
Results are qualitatively similar in both cases suggesting that the preparation of LDA does not 
alter our conclusions. 
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